GEOG 331: Environmental Data Science, Colgate University

Instructions

There are 13 questions to this activity. Save your answers in word document that you will hand in on Moodle using a .pdf extension. Keep your script file in your GitHub folder and make sure that all changes are pushed to GitHub. Code for all questions should be clearly commented with the question number. You will include a link to your script file as a part of your final question in this activity. The assignment is worth 50 points.

Learning objectives

1. Evaluate glacial recession and vegetation change in Glacier National Park

2. Learn to work with vector and raster data in R

3. Work with projections in R

4. Conduct an analysis of spatial data

5. Make a nice map

Introduction to the problem: glacial retreat and vegetation change

As air temperatures rising under global climate change, areas covered in snow and ice are melting. Glaciers are masses of ice that are currently melting as air temperatures increase. Glaciers form when more snow accumulates than melts or evaporates and the snow becomes compacted over time. Glaciers are defined as being persistent over time (don’t disappear seasonally or even every few years) and must be at least 0.1 km squared in area, and often take thousands of years to form. Glaciers often form in mountainous regions. Glacier National Park is located in Montana and known for sweeping views of glacier covered mountains. The glaciers in Glacier National Park formed at least 7,000 years ago. However, the glaciers in Glacier National Park began noticeably shrinking over the last century. There is concern about impacts on tourism, ecosystems, and hydrology as the glaciers disappear. As shown below, the receding glaciers are changing the landscape.

Source: National Park Service


Glacial retreat will drastically affect the regional hydrology, vegetation, and the energy balance. For example, a surface with snow or a glacier reflects more incoming solar radiation off of the earth’s surface compared to bare rock and vegetation nearby. The loss of glaciers means that the land surface will experience even greater warming as more solar energy is absorbed. As glaciers retreat and overall snow and ice cover decreases in mountainous regions, vegetation can grow in newly formed bare surfaces. You can see in the picture above, there is grass growing in areas that were once covered in snow and ice. This growth of vegetation can drive additional change on the land surface and alter water, energy, and carbon cycling. A study by Carlson et. al. in 2017 demonstrated that less snow cover during the year drives increased vegetation in the French Alps. Another recent study by Karen Anderson et. al. in 2019 found that the Himalayas are experiencing increased vegetation growth as glaciers and permanent snow cover retreats. Read through these papers posted in the data folder to get an idea of the type of analysis and considerations for this type of research.

In this exercise, you’ll use data from the National Park Service that shows the footprint of 39 glaciers in 1966, 1998, 2005, and 2015 to show the change in glacier extent over recent decades. We’ll also look at a satellite derived measure of vegetation called the Normalized Difference Vegetation Index (NDVI). NDVI Values ranging from 0 to 1 indicate the relative amount of photosynthetic material on the ground. Values close to one indicate a higher level of vegetation and zero indicates no vegetation. The overall value will depend on characteristics of vegetation, and areas completely covered in vegetation do not necessarily correspond to an NDVI value of one. NDVI varies seasonally and typically peaks in the summer for many cold winter ecosystems. The maximum NDVI value indicates the highest amount of vegetation that year. Tracking changes in NDVI throughout time can indicate if there is increasing or decreasing amounts of vegetation.

Set up spatial packages

You will use a series of spatial packages. Install the following packages:

#install.packages(c("raster","sp","rgdal","rgeos","plyr"))
You’ll notice there is a package plyr on the list. It is an older version of the package dplyr. We can’t use dplyr here because it overwrites too many spatial functions. If you were to load dplyr, you would get the following messages:


That is a lot of useful spatial functions that are overwritten by dplyr functions. You definitely want to avoid this when working with spatial data. Installation may take a few minutes. Next load all the packages to your environment.

library(raster)
library(sp)
library(rgdal)
library(rgeos)
library(plyr)

Reading in vector data

You will start by reading in the vector data of the glacier outlines. These glaciers have been digitized by the National Park Service. In shapefiles, all metadata is contained in the .xml file. You can open it in a web browser to find more information. The GNP_glaciers file contains all four shapefiles for each year. Read them in using the readOGR function.

#read in shapefiles
#readOGR in rgdal does this
g1966 <- readOGR("data\\GNPglaciers\\GNPglaciers_1966.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_1966.shp", layer: "GNPglaciers_1966"
## with 39 features
## It has 13 fields
## Integer64 fields read as strings:  OBJECTID
g1998 <- readOGR("data\\GNPglaciers\\GNPglaciers_1998.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_1998.shp", layer: "GNPglaciers_1998"
## with 39 features
## It has 13 fields
## Integer64 fields read as strings:  OBJECTID
g2005 <- readOGR("data\\GNPglaciers\\GNPglaciers_2005.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_2005.shp", layer: "GNPglaciers_2005"
## with 39 features
## It has 13 fields
## Integer64 fields read as strings:  OBJECTID
g2015 <- readOGR("data\\GNPglaciers\\GNPglaciers_2015.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_2015.shp", layer: "GNPglaciers_2015"
## with 39 features
## It has 21 fields
## Integer64 fields read as strings:  OBJECTID Recno

Note there are 39 features (glaciers) for each shapefile and a table with information about each glacier that has 13 columns of data. Let’s investigate what the format of this data.

str(g2015)
#data stores all accompanying info/measurements for each spatial object
head(g2015@data)
##   OBJECTID Recno Year Src1_0822 Src1_nadir Src3_nadir Src2_nadir   Area2015
## 0        1  1055 2015 secondary    11.7240     0.0000     4.7346  736669.75
## 1        2  1202 2015   primary    27.0008     0.0000     0.0000  511589.79
## 2        3  1120 2015 secondary    11.7240     0.0000     4.7346   75562.60
## 3        4  1430 2015 secondary     0.0000    16.8225     0.0000 1498505.92
## 4        5  1041 2015   primary    11.7240     0.0000    26.2766   35298.01
## 5        6  1132 2015 secondary    11.7240     0.0000     4.7346  224773.89
##   Shape_leng  X_COORD Y_COORD                                            SOURCE
## 0   7741.335 268429.9 5425167 WorldView-01 Satellite imagery from Digital Globe
## 1   5184.650 295750.8 5413701 WorldView-01 Satellite imagery from Digital Globe
## 2   1255.828 268868.3 5421731 WorldView-01 Satellite imagery from Digital Globe
## 3  15290.128 303278.4 5385710 WorldView-01 Satellite imagery from Digital Globe
## 4   2259.001 273823.9 5427266 WorldView-01 Satellite imagery from Digital Globe
## 5   4090.754 276459.4 5419663 WorldView-01 Satellite imagery from Digital Globe
##             CLASSIFICA             OWNERSHIP SrcOth_nad SrcOth_dat
## 0 main body of glacier Glacier National Park     0.0000 9999/01/01
## 1 main body of glacier Glacier National Park    18.6183 2014/10/19
## 2 main body of glacier Glacier National Park     0.0000 9999/01/01
## 3 main body of glacier Glacier National Park     0.0000 9999/01/01
## 4 main body of glacier Glacier National Park     0.0000 9999/01/01
## 5 main body of glacier Glacier National Park     0.0000 9999/01/01
##                                                               COMMENT Src2_0912
## 0            used both 20150822 and 20150912 images for full coverage   primary
## 1  20150822 only 2015 image, but high nadir, also used 20141019 image       n/a
## 2             20150822 used where shading was heavy in 20150912 image   primary
## 3   used color 20150822 as secondary, some offset from 20150925 image       n/a
## 4 used image from 20150822 since later image had high off nadir angle secondary
## 5                                                          no comment   primary
##   Src3_0925          GLACNAME SOURCE_SCA
## 0       n/a   Agassiz Glacier    1:12000
## 1       n/a     Ahern Glacier    1:12000
## 2       n/a      Baby Glacier    1:12000
## 3   primary Blackfoot Glacier    1:12000
## 4       n/a   Boulder Glacier    1:12000
## 5       n/a    Carter Glacier    1:12000

This table allows for data describing a spatial object to be stored. You can work with this table just like any other table, but you will be able to display objects in it spatially.

#polygons stores the coordinates for drawing the polygons
g2015@polygons[[1]]

You’ll notice there can be multiple polygons within a polygon. This is because there may be multiple, unconnected shapes that are considered to have the same characteristics and share the same data in the table. You’ll also see all the items needed for drawing a polygon including the coordinates, any holes to put in the polygon, and coordinates for the center. For a vector object, you can find the projection info in an object called proj4string

g1966@proj4string
## CRS arguments:
##  +proj=utm +zone=12 +datum=NAD83 +units=m +no_defs


Now let’s make a map of the glaciers and view them by name.The spplot function allows you to map vector data and show different colors for a data value. Let’s look at the 1966 glaciers by glacier name.

spplot(g1966, "GLACNAME")

Before you get started, there’s one other issue. The glacier names in 2015 don’t match the rest of the shapefiles. Compare 2015 to 1966 for reference:

#check glacier names
g1966@data$GLACNAME
##  [1] Agassiz Glacier         Ahern Glacier           Baby Glacier           
##  [4] Blackfoot Glacier       Boulder Glacier         Carter Glacier         
##  [7] Chaney Glacier          Dixon Glacier           Gem Glacier            
## [10] Grant Glacier           Grinnell Glacier        Harris Glacier         
## [13] Harrison Glacier        Herbst Glacier          Hudson Glacier         
## [16] Ipasha Glacier          Jackson Glacier         Kintla Glacier         
## [19] Logan Glacier           Lupfer Glacier          Miche Wabun Glacier    
## [22] N. Swiftcurrent Glacier Old Sun Glacier         Piegan Glacier         
## [25] Pumpelly Glacier        Rainbow Glacier         Red Eagle Glacier      
## [28] Salamander Glacier      Sexton Glacier          Shepard Glacier        
## [31] Siyeh Glacier           Sperry Glacier          Stanton Glacier        
## [34] Swiftcurrent Glacier    Thunderbird Glacier     Two Ocean Glacier      
## [37] Vulture Glacier         Weasel Collar Glacier   Whitecrow Glacier      
## 39 Levels: Agassiz Glacier Ahern Glacier Baby Glacier ... Whitecrow Glacier
g2015@data$GLACNAME
##  [1] Agassiz Glacier            Ahern Glacier             
##  [3] Baby Glacier               Blackfoot Glacier         
##  [5] Boulder Glacier            Carter Glacier            
##  [7] Chaney Glacier             Dixon Glacier             
##  [9] Gem Glacier                Grant Glacier             
## [11] Grinnell Glacier           Harris Glacier            
## [13] Harrison Glacier           Herbst Glacier            
## [15] Hudson Glacier             Ipasha Glacier            
## [17] Jackson Glacier            Kintla Glacier            
## [19] Logan Glacier              Lupfer Glacier            
## [21] Miche Wabun                North Swiftcurrent Glacier
## [23] Old Sun Glacier            Piegan Glacier            
## [25] Pumpelly Glacier           Rainbow Glacier           
## [27] Red Eagle Glacier          Salamander Glacier        
## [29] Sexton Glacier             Shepard Glacier           
## [31] Siyeh Glacier              Sperry Glacier            
## [33] Stanton Glacier            Swiftcurrent Glacier      
## [35] Thunderbird Glacier        Two Ocean Glacier         
## [37] Vulture Glacier            Weasel Collar Glacier     
## [39] Whitecrow Glacier         
## 39 Levels: Agassiz Glacier Ahern Glacier Baby Glacier ... Whitecrow Glacier

You’ll notice Miche Wabun is missing Glacier in its name and North Swiftcurrent Glacier should be abbreviated as N.. Fix that so you can use the glacier names as an indexing variable across years.

#fix glacier name so that it is consistent with the entire time period
g2015@data$GLACNAME <- ifelse(g2015@data$GLACNAME == "North Swiftcurrent Glacier",
                              "N. Swiftcurrent Glacier",
                                          ifelse(   g2015@data$GLACNAME ==  "Miche Wabun", 
                                               "Miche Wabun Glacier",
                                                as.character(g2015@data$GLACNAME)))

Working with raster data

Mapping satellite imagery

It’s hard to get a reference for the glaciers just from the polygon file. Let’s read in an image collected by the Landsat satellite on September 5, 2014. This will provide some reference for the location. The landsat sensors measure the light that is reflected up to the top of the atmosphere. If we overlay the red, green, and blue visible bands of light, it will look just like a photograph of the earth surface. Here you will work with the raster package to read in these .tif files. These files are different from regular photograph tif files since they are actually geotiff files that contain spatial reference information.

#read in rgb imagery from landsat
redL <- raster("data\\glacier_09_05_14\\l08_red.tif")
greenL <- raster("data\\glacier_09_05_14\\l08_green.tif")
blueL <- raster("data\\glacier_09_05_14\\l08_blue.tif")

You can check the coordinate system information in the raster files by using the following code:

#check coordinate system
redL@crs
## CRS arguments:
##  +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs

You’ll notice the coordinate system is the same as the vector data above. This means that you can plot them together without having to reproject. It will be easier to plot these rasters together if we turn them into a brick. A brick is a series of raster files with the same extent and resolution. This object allows us to readily overlay the objects.

#make a brick that stacks all layers
rgbL <- brick(redL, greenL, blueL)

Now that the raster files are grouped together, make a map with both the polygons. You’ll use plotRGB to make a color image of the raster. Using the plot function adds the polygons to the map. Notice the units of the x and y axis here. When you plot polygons, you can specify a single color for a polygon or if you create a vector of colors that is the same length of the total number of polygon objects (39 here).

#plot with color
#show axes for reference
#add contrast to the imagery to see it better
par(mai=c(1,1,1,1))
plotRGB(rgbL, stretch="lin", axes=TRUE)
#add polygons to plot
plot(g1966, col="tan3", border=NA, add=TRUE)
plot(g1998, col="royalblue3", add=TRUE, border=NA)
plot(g2005, col="darkgoldenrod4", add=TRUE, border=NA)
plot(g2015, col="tomato3", add=TRUE, border=NA)

It is a little hard to see the details at the scale of the entire park. You can use the ext argument to zoom in on a specific area. ext relies on the following order: xmin, xmax, ymin, ymax.Subset the plot below to zoom in on closer on a few glaciers:

plotRGB(rgbL, ext=c(289995,310000,5371253,5400000), stretch="lin")
plot(g1966, col="palegreen2", border=NA, add=TRUE)
plot(g1998, col="royalblue3", add=TRUE, border=NA)
plot(g2005, col="darkgoldenrod4", add=TRUE, border=NA)
plot(g2015, col="tomato3", add=TRUE, border=NA)

Working with raster data

The NDVI folder contains NDVI data collected from the MODIS sensor from 2003-2016. The NDVI data is the maximum NDVI that is observed for each year and is a derived data product generated by the USGS eModis Phenology program. The maximum NDVI will characterize the peak level of vegetation during the growing season. Let’s read in all of these files. We’ll read them in as a list and each item in our list will be one year of raster data.

#set up years to read in
ndviYear <- seq(2003,2016)

#read all files into a list
NDVIraster <- list() 
for(i in 1:length(ndviYear)){
    NDVIraster[[i]] <- raster(paste0("data\\NDVI\\NDVI_",ndviYear[i],".tif"))

}
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition

Use the str function to take a look at what is in a single raster function. Look at the first year (2003):

str(NDVIraster[[1]])

Take a closer look at the coordinate system:

#get projection
NDVIraster[[1]]@crs
## CRS arguments:
##  +proj=laea +lat_0=48 +lon_0=-113 +x_0=0 +y_0=0 +ellps=sphere +units=m
## +no_defs
This is a Lambert Azimuthal Equal Area projection that has an orgin (centers the zero point of the map at that latitude and longitude) of 48 degrees and -113 degrees. x_0 and y_0 indicate if there is a false easting or northing set (moves zero point in coordinate system, typically to keep all numbers positive). +a and +b describe the radius of the ellipsoid axes used for the datum. Units indicate the units the project is in. Many projections are in meters, but some US projections will be in feet. Note that any additional latitude arguments (lat_1, lat_2) would be associated with a standard parallel (lines of high accuracy in some projections). There are none here. This is going to be a good projection to work in since it is an equal area projection.


Let’s take a closer look at the NDVI data by making a plot of it:

plot(NDVIraster[[1]])

You’ll notice there are a few values below zero for NDVI. Values below zero represent water and other non-vegetative features. For this analysis we will ignore them. If you were using this data for scientific research, you might want to better filter the data.


## Vector data analysis: glacier retreat The first step to analyzing glacial retreat is to make sure the glaciers are in the right projection. The NDVI raster data is in a good projected coordinate system (PCS), and we can project the glaciers into the same PCS. The projection tool will be in the sp package for vector data (note rasters have their own projection tool in the raster package). You can just directly reference the NDVI raster projection information.

#reproject the glaciers
#use the NDVI projection
#spTransform(file to project, new coordinate system)
g1966p <- spTransform(g1966,NDVIraster[[1]]@crs)
g1998p <- spTransform(g1998,NDVIraster[[1]]@crs)
g2005p <- spTransform(g2005,NDVIraster[[1]]@crs)
g2015p <- spTransform(g2015,NDVIraster[[1]]@crs)


Note: that sometimes resizing plots with multiple sources of spatial data can have issues resizing all data layers. If it looks like the extent doesn’t match after resizing, rerun your code again after resizing.

Next let’s analyze the area of the glaciers throughout each year. There is an area column in the glacier data tables, but the metadata does not indicate how the area was measured. We will calculate the area directly for each polygon so that we know that we are consistently measuring the area of each polygon. The area function will return a vector of areas in our coordinate system units for each of the polygons.

#calculate area for all polygons
#add directly into data table for each shapefile
g1966p@data$a1966m.sq <- area(g1966p)
g1998p@data$a1998m.sq <- area(g1998p)
g2005p@data$a2005m.sq <- area(g2005p)
g2015p@data$a2015m.sq <- area(g2015p)

It will be helpful to join all of this data together into a table not associated with the shapefile so that you can analyze it. Here you will use the join from plyr.

gAllp1 <- join(g1966p@data,g1998p@data, by="GLACNAME", type="full")
gAllp2 <- join(gAllp1,g2005p@data, by="GLACNAME", type="full")
gAll <- join(gAllp2,g2015p@data, by="GLACNAME", type="full")

Let’s make a plot of the area for each glacier using this table:

plot(c(1966,1998,2005,2015), 
        c(gAll$a1966m.sq[1],gAll$a1998m.sq[1], gAll$a2005m.sq[1],gAll$a2015m.sq[1]),
        type="b", 
        pch=19, col=rgb(0.5,0.5,0.5,0.5), xlim= c(1965,2016),
        ylim=c(0,2000000),
        ylab="Area of glacier (meters squared)",
        xlab="Year")
        
for(i in 2:39){
points(c(1966,1998,2005,2015), 
        c(gAll$a1966m.sq[i],gAll$a1998m.sq[i], gAll$a2005m.sq[i],gAll$a2015m.sq[i]),
        type="b", 
        pch=19, col=rgb(0.5,0.5,0.5,0.5))

}

This is a good first look, but it would be more helpful to look at the percent change since small glaciers overall have less area to lose.


It might help to visualize the glacial loss more directly. Let’s make a polygon that shows the difference in glaciers between 2015 and 1966. Here, you can use the many vector operations available in rgeos. Since we want to view the area where 1966 and 2015 don’t overlap, the gDifference function will remove any areas that overlap.

diffPoly <- gDifference(g1966p, g2015p, checkValidity = 2L)
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -79835.564606180007 107115.29284187
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -77011.733007000003 99328.148267840006
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -59187.508340740002 93811.335311250004
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -62514.363778120001 96430.893114039995
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -75367.583780639994 92136.2396056
## g2015p is invalid
## Attempting to make g2015p valid by zero-width buffering
plot(diffPoly)

What does the checkValidity argument in the gDifference function do? Why do you think we used that?

Let’s look at the difference with more context now. Let’s plot the glacier retreat with the map of NDVI showing all areas of the glacial loss in black:

#plot with NDVI
plot(NDVIraster[[13]], axes=FALSE, box=FALSE)
plot(diffPoly,col="black", border=NA,add=TRUE)

This difference polygon is nice for plotting, but if you run str you will notice it lost the original geometry and linked data. You may want to work with the orignial data if you need to maintain glacial names and other information.


Raster data analysis: does more vegetation grow with glacial retreat?

Let’s take a look at patterns in maximum NDVI from 2003-2016 to see if there is an increase in vegetation throughout the time period. Let’s start by examining how NDVI might be changing within the area that glaciers are retreating. We can get all of the values that are within a polygon using the extract function. Let’s extract each year’s NDVI in the area that experienced retreat from 1966 to 2015 and take the average. We’ll have to use a for loop to work with all values in the dataset.

#extract NDVI values
NDVIdiff <- list()
meanDiff <- numeric(0)
#loop through all NDVI years
for(i in 1:length(ndviYear)){
  #get raster values in the difference polygon
    NDVIdiff[[i]] <- extract(NDVIraster[[i]],diffPoly)[[1]]
    #calculate the mean of the NDVI values
    meanDiff[i] <- mean(NDVIdiff[[i]], na.rm=TRUE)
}

Let’s take a look at the data:

plot(ndviYear, meanDiff, type="b",
    xlab= "Year",
    ylab="Average NDVI (unitless)",
    pch=19)

Although 2015 had very high NDVI, it doesn’t appear that there is much of a trend.This is not surprising because it may take many years for vegetation to move in and fill in a newly bare glacial area. There may also be prolonged snow cover in many of these areas even if there is not a glacier there. Let’s see if vegetation is increasing in areas surrounding the glaciers. First it will be helpful to know how much NDVI is changing in each pixel. We can fit a slope using the lm function to each cell in the pixel. We can use the calc function which applies a function to every pixel in the raster. Since you will want to apply the function to every cell in each year and for every yearly raster, you will want to be sure that R knows the raster is a stack (similar to brick) of rasters.

#designate that NDVIraster list is a stack
NDVIstack <- stack(NDVIraster)
#set up lm function to apply to every cell
#where x is the value of a cell
#need to first skip NA values (like lakes)
#if NA is missing in first raster, it is missing in all
#so we can tell R to assign an NA rather than fitting the function
timeT <- ndviYear
fun <- function(x) {
    if(is.na(x[1])){
        NA}else{
        #fit a regression and extract a slope
            lm(x ~ timeT)$coefficients[2] }}
#apply the slope function to the rasters
NDVIfit <- calc(NDVIstack,fun)
#plot the change in NDVI
plot(NDVIfit, axes=FALSE)

The raster NDVI fit now contains the slope for the change in NDVI per year.


Let’s examine how this change might be occurring around the glaciers. We can use the gbuffer function to look at the area that is 500 meters out from the largest glacier extent to make sure we are outside of the glacier. Let’s examine if there are any individual differences in each glacier.

#buffer glaciers
glacier500m <- gBuffer(g1966p,#data to buffer
                    byid=TRUE,#keeps original shape id 
                    width=500)#width in coordinate system units

We can use zonal statistics to get a statistic applied across an index, but first we need to convert our vector data to raster. The zonal function only works with two rasters.

#convert to a raster
buffRaster <- rasterize(glacier500m,#vector to convert to raster
                    NDVIraster[[1]], #raster to match cells and extent
                    field=glacier500m@data$GLACNAME, #field to convert to raster data
                    background=0)#background value for missing data
plot(buffRaster)

We also need to remove the actual glacier from our statistics. Let’s rasterize the largest glacial extent. Since the IDs will be the same, we can subtract the two rasters to end up with an id for only in areas with the buffer.

#rasterize gralciers
glacRaster <- rasterize(g1966p, NDVIraster[[1]], field=g1966p@data$GLACNAME, background=0)
#subtract buffer from original glacier
glacZones <- buffRaster - glacRaster
plot(glacZones)


Now you can get the statistics for the rate of vegetation change in the area around the rasterized buffer using zonal statistics:

meanChange <- zonal(NDVIfit, #NDVI function to summarize
                glacZones,#raster with zones
                "mean")#function to apply
head(meanChange)
##      zone          mean
## [1,]    0  0.0011640432
## [2,]    1  0.0008101119
## [3,]    2  0.0001584056
## [4,]    3  0.0018559251
## [5,]    4 -0.0014669432
## [6,]    5  0.0027154395